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Abstract. The control of quantum dynamics via specially tailored laser pulses is 
a long-standing goal in physics and chemistry. Partly, this dream has come true, as 
sophisticated pulse shaping experiments allow to coherently control product ratios of 
chemical reactions. The theoretical design of the laser pulse to transfer an initial state 
to a given final state can be achieved with the help of quantum optimal control theory 
(QOCT). This tutorial provides an introduction to QOCT. It shows how the control 
<**] ■ equations defining such an optimal pulse follow from the variation of a properly defined 

functional. We explain the most successful schemes to solve these control equations and 
^-i show how to incorporate additional constraints in the pulse design. The algorithms are 

then applied to simple quantum systems and the obtained pulses are analyzed. Besides 
the traditional final-time control methods, the tutorial also presents an algorithm and 

■ an example to handle time-dependent control targets. 
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■ 1. Introduction 

H . 

Since the first realization of the laser by T.H. Maiman in 1960, physicists and chemists 
have had the vision to coherently control quantum systems using laser fields. For 
example, laser pulses may be applied to create and break a particular bond in a molecule, 
to control charge transfer within molecules, or to optimize high harmonic generation. 
The first approach to break a certain bond in a molecule using a laser tuned on resonance 
and initiating a resonance catastrophe failed [I]. The molecule internally converted the 
energy too quickly, so that the specific bond did not break but instead the whole molecule 
was "heated" [2]. To overcome this so-called internal vibrational relaxation (IVR) a 
smarter excitation strategy and further technological improvements were necessary. 

With the advent of femtosecond laser pulses in the 1980's and a sophisticated pulse- 
shaping technology [3] the goal of controlling complex chemical reactions with coherent 
light was finally achieved: For example, in 1998 Assion et al. pi] showed that the product 
ratio CpFeCOCl + /FeCl + of the organo-metallic compound (CpFe(C0)2Cl) can be either 
maximized or minimized by a specially tailored light pulse; or in 2001, Levis et al. [5] 
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demonstrated a rearrangement of molecular fragments. In both of these experiments 
adaptive laser pulse shaping techniques [3l|6] were applied, i.e., a computer analyzes the 
outcome of the experiment and modifies the laser pulse shape to optimize the yield of 
a predefined reaction product. This process is repeated until the optimal laser pulse is 



found (see Appendix A.l ). The number of experiments based on this so-called closed- 



learning- loop (CLL) is growing constantly, see for example Refs. [7] l8l |9| [TUl [TT| [12l [13] . 
Recently, the pulse shaping techniques have been extended to allow also for polarization 
shaping [HI [151 [16j Q2], i.e., experimentalists can independently shape polarization, 
amplitude, and phase. 

In addition to further technological advances, it is of utmost importance to have 
powerful theoretical methods available. The questions that theory has to answer can be 
divided into two classes: The first class is that of controllability [18j, or in other words: 
Given a certain quantum system (e.g. a molecule) can the control target (a certain 
reaction product) be reached at all with the given controller (e.g. a laser)? The second 
class concerns the problem of finding the best way to achieve a given control objective, 
e.g., calculating the optimal laser pulse for breaking a particular bond in a molecule. 
Such theoretical predictions are very important to gain insight into the complexity of 
the control process, to determine the experimental parameters, to transfer laser pulse 
designs to the laboratory, or to compare the optimized pulse from an experiment with 
the calculated one [19J. 

Several theoretical approaches have been developed to optimize laser pulses, ranging 
from brute- force optimization of a few pulse parameters [20j, pulse-timing control 
[2U [221 [23] , Brumer-Shapiro coherent control [21], stimulated- Raman- Adiabatic-Passage 
(STIRAP) [25l EEl [27] , to genetic algorithms [28J. The most powerful approach, in our 
opinion, is optimal control theory (OCT) which is commonly applied in engineering, for 
example to design trajectories for satellites and space probes. The application of optimal 
control theory to quantum mechanics started in the late 1980 's [291 E0] and shows 
continuous advances until today. Among the most important developments were the 
introduction of rapidly converging iteration schemes [311 E21 [33], El] , the generalization 
to include dissipation (Liouville space) [351 EHl [37] , and to account for multiple control 
objectives [38] . 

This tutorial will focus on the theoretical aspects of quantum optimal control theory 
(QOCT) and tries to explain the beginning Phd student how to calculate optimal laser 
pulses. The first step is to work through the basic theory which is presented in section EJ 
In the following sections we will always refer to parts of the presented theory therein. 
The application of the derived method is then explored with the help of a two-level 



system in section 12.71 The basics of the two- level system are reviewed in Appendix A. 2 



As a second example we consider the control of a ID asymmetric double well model 
in section 12.81 which will motivate the need for further constraints on the optimal laser 
field. Theory and algorithms for optimizing laser pulses with additional constraints are 
explained in section [31 Two examples for the asymmetric double well can be found in 
section 13.41 and section 13.51 The last part of this tutorial in section 0] will focus on time- 
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dependent control targets. A brief summary and outlook of this article can be found in 
section O 

We would like to conclude this introduction by emphasizing that our intention 
is to provide only a brief overview, but a detailed derivation of the algorithms, and 
simple examples that can be reworked by the reader (especially those of the two-level 
system). The tutorial cannot cover all topics in quantum control and excludes the 
following imported topics: Closed- loop control (see Refs. [39|), time optimal control 
(see Refs. pH SH H2] ) , and dissipative systems (see Refs. [351 ES M, 1131 Hi])- For an 
excellent review on the experimental aspects of quantum control we like to refer the 
reader to Ref. [35] . 

2. Theory and Algorithms 

In this section we sketch the basics of optimal control theory applied to quantum 
mechanics. 

Let us consider an electron in an external potential V(r) under the influence of a 
laser field propagating in z-direction. Given an initial state \l/(r, 0) = 0(r), the time 
evolution of the electron is described by the time-dependent Schrodinger equation with 
the laser field modelled in the dipole approximation (length gauge) 

ij^(r,t) = jftt(r,t), (1) 

H = H -fie(t), (2) 

H =f+V, (3) 

(atomic units are used throughout: h — m — e — 1). Here, fi = (£l x , jl y ) is the dipole 
operator, and e(t) = (e x (t) , e y (t)) is the time-dependent electric field. The kinetic en- 
ergy operator is T = — 

2.1. Controllability 

Before trying to find an optimal control for a given target and quantum system we 
raise the question if the given control target can be reached at all with the given 
controller, here: the controller is the laser field Hi = —fie(t). In the following we want 
to summarize some of the rigorous results on controllability that exist in the literature. 

The most powerful and easy-to-use statements are available for A^-level systems 
[T8l HHl H?]. We start with a definition of the term complete controllability and then 
discuss the results from Ref. [16] . 

Definition (Schirmer et al. fj^j )-' A quantum system 



H =H + Hj, (4) 

N M 

H = J2 £ n\n)(n\, H I = J2frn{t)H m , (5) 

n=l m=l 



Quantum Optimal Control Theory 



4 



is completely controllable if every unitary operator U is accessible from the identity 
operator J via a path j(t) = U(t,t Q ) that satisfies 

id t U(t,t ) = (H + Hi)u(t,t ). (6) 



Theorem (Ramakrishna et al. [^6]): A necessary and sufficient condition for 
complete controllability of a quantum system defined by equations (HI) and ([5]) is that 
the Lie algebra L Q has dimension N 2 . 

We can therefore check the controllability of an iV-level system by constructing 
its Lie algebra L which is generated by Ho, ... , Hm and then calculate the rank of 
the algebra. An algorithm for this task has been suggested in Ref . [47] . This scheme is 



demonstrated for the two- level-system in Appendix A .4 A further statement in Ref. [46 



guarantees that complete controllability can also be achieved under external constraints 
on the strength of the controller. 

The extension of such controllability theorems to an infinite-dimensional Hilbert 
space and the inclusion of unbound operators like r or V r turns out to be non-trivial as 
shown in Ref. [18] . The conditions of controllability are only valid for quantum systems 
with a non-degenerate and discrete spectrum and do not include external constraints 
on the strength of the control functions. 

2.2. Derivation of the control equations 

Let us consider the following quantum mechanical control problem: 

Our goal is to find a laser pulse e(t) which drives a quantum system from its initial 
state ^(O) to a state *S?(T) in such a way that the expectation value of an operator O is 
maximized at the end of the laser interaction: 



maxJ! with Ji[tf] = (tf(T)|0|tf(T)). (7) 

e(t) \ / 

At this point we keep the operator O as general as possible. The only restriction on 
O is that it has to be a Hermitian operator. We will discuss examples for the operator 
O in section [231 In addition to the maximization of we require that the fluence 

of the laser field is as small as possible which is cast in the following mathematical form: 

f T 

M e }= - z2 / dt a 3 e j 2 ^' j = x ,y- ( 8 ) 

where e x (t) and e y {t) are the components of the laser field perpendicular to the 
propagation direction. The positive constants (Xj play the role of penalty factors: The 
higher the laser fluence the more negative the expression, the smaller the sum Ji + J 2 . As 
pointed out in Ref. [IE] , the penalty factor can be extended to a time-dependent function 
atj(t) to enforce a given time- dependent shape of the laser pulse, e.g. a Gaussian or 
sinusoidal envelope. The constraint that the electronic wave function has to satisfy the 
time-dependent Schrodinger equation (TDSE) is expressed by 

J 3 [e,*, X ]= ~ 2hn J dt (x(t)|(i$-tf(t))|*(t)), (9) 
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where we have introduced the Lagrange multiplier x(t)- Since we require the TDSE 
to be fullfilled by the complex conjugate of the wave- function as well, we obtain the 
imaginary part Im of the functional. Note, that we choose the imaginary part to be 
consistent with the literature, for instance, Ref. [301 H6| [32]. 
The Lagrange functional has the form 

J[ X , , e] = + J 2 [e] + J 3 [x, , e]. (10) 

We refer to this functional as the standard optimal control problem and start the 
discussion of all extensions considered in this work from this standard form. 



2. 3. Variation of J 



To find the optimal laser field from the functional in equation (flOj) we perform a total 
variation. Since the variables x an d e are linearly independent we can write: 



5J 



drldr 
o J U*(r,r 



3 J 6 J 1 

5V(t,t) + - f _ _^ 5 X (t,t) \ + J2 / dT 

) k=x,y J ° 



5J 



5eu(t) 



= 5^J + 5 X J + 2J ^e fc ^> 
k=x,y 

where we have omitted the derivations with respect to the complex conjugate of the 
wave function ty*(t) and the Lagrange multiplier x*(0> since the functional derivative 
will result in the complex conjugate equations for these variations. 

Since we are looking for a maximum of J, the necessary condition is 

5J = 

=^<S v J = 0, <LJ = 0, 5 efc ^ = 0. (12) 



-Se k (r) 



(11) 



2.3.1. Variation with respect to the wave function $ First, let us consider the functional 
derivative of J with respect to ^: 

-J£—=d**(r',T)8(T-T), 
oW(r , r) 

— = 

sj 3 

For the last functional derivative we have used the following partial integration: 



- i (id T + H(t) X *(r', r) - [ X *(r', t)5(t - r)} 



dt ( X (t)\idt-H(t)Mt) 



i(v(/)|*f-/-)) T -i I dt (d tX (t)m)) 
Jo 







dt (H{t)x{t)Mt) 



dt ((id t -H(t))x(t)mt))- 



(13) 



Quantum Optimal Control Theory 
We find for the variation with respect to 

6*J= (*(T)\6\6#(T))+i J dT((id T -H{r))x(T)\S*(T) 

-( X (T)|^(T)) + (x(0)|^(0)). (14) 

S V ' 

=0 

The variation of <5\I/ (0) vanishes because we have a fixed initial condition, \l/(0) = 0j. 

2.3.2. Variation with respect to the Lagrange multiplier x Now we do the same steps 
for x(t)- 

5Ji 5J 2 



Sx(r',r) 5x(r',r 
SJ 3 



0. 



i(i0 T + #(r)) r(r',r). 



Sx{r',r) 

In contrast to the variation with respect to \l/ we do not have boundary terms here. The 
variation of J with respect to x yields 

S X J= -i J dr^(idT-H(T)^{T)\6x(r)), k = x,y. (15) 

2.3.3. Variation with respect to the field The functional derivative with respect to €k(t) 
is 



5e k {r) 



0. 

-2a fc e fc (r), (16) 
-2Im(x(r)|/i fc |*(r)), k = x,y. (17) 



<$e fc (r) 

Hence, the variation with respect to €k(t) yields 



5eJ= I dr {-21m{ X {r)\fl k \^{T))-2a k e k {T)}5e k {T). 
Jo 



;is) 



2.4- Control equations 

Setting each of the variations independently to zero results in the desired control 
equations. From 5 ek J = [equation ([IB]) ] we find 

a k e k (t) = -lm( X (t)\ix k \^(t)), k = x,y. (19) 

The laser field e(t) is calculated from the wave function \l/(t) and the Lagrange multiplier 
x{t) at the same point in time. The variation 5 X J in equation f|T5]) yields a time- 
dependent Schrodinger equation for \l/(t) with a fixed initial state (pi, 

= (id* - *(r, t), *(r, 0) = &(r). (20) 

Note that this equation also depends on the laser field e(t) via the Hamiltonian. 
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The variation with respect to the wave function 8^ J in equation (1141) results in 

(id t - X (r, t) = i (x(r, t) - 6*(r, t)) 5(t - T). (21) 

If we require the Lagrange multiplier x(t) to be continuous at t — T, we can solve the 
following two equations instead of equation (j2~Tj) : 

'idt-H(tj) X (T,t) = 0, (22) 

X (r,T) =0*(r,T). (23) 

To show this we integrate over equation ( 1211 : 

(id t -A(t)) X (r,t) 

5(t-T). (24) 

The left-hand side of equation (1241) vanishes because the integrand is a continuous 
function. It follows that also the right-hand side must vanish, which implies equation 
(1231) . From equations ( 1231) and ( 12TT) then follows equation ( 1221) . Hence, the Lagrange 
multiplier x(t) satisfies a time- dependent Schrodinger equation with an initial condition 
at t = T. The set of equations that we need to solve is now complete: ( fl9l) . ( 1201) . ( 1221) . 
and ([23]). 

To find an optimal field e(i) from these equations we use an iterative algorithm 
which will be discussed in the section 12.61 





pT+K f 


lim 


/ dt 




It-k - 




rT+K 


lim 


/ dti 




IT-K 



2.5. Target operators 

So far we have shown how we can optimize a laser field and which equations have to be 
solved to achieve this goal. Before we discuss details on how the equations are solved in 
practice, we want to discuss different examples for the target operator O. 

2.5.1. Projection operator Choosing a projection operator O = \<l>f)(<f>f\, the 
maximization of J\ corresponds to maximizing the overlap of the propagated wave 
function ty(T) with <ftf, i.e., 

Ji = (*(T)|^> (0/|*(T)) = ((^(T)!^)! 2 . (25) 

Using a projection operator as target in the optimal control algorithm therefore allows 
one to find an optimized pulse which drives the system from the initial state ^(0) to 
the desired target state 0/ up to a global phase factor e 17 . In other words, the target 
functional J\ is invariant under the transformation 

4>t-+&4>f. (26) 
It is possible to fix this phase if we replace equation (1231) by the following functional 

He(*(T)|^>, (27) 
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which can be derived in the following way: 

min ||^(T) - <P f \\ 2 = min {(tf (T)|tf (T)) + (0/(0/) - 2 Re (V(T)\4> f )} (28) 

Assuming normalized states the minimization corresponds to the maximization of 
equation (1271) . 

In all cases considered in this work, 0/ is chosen to be an excited state of the 
quantum system. But it is also possible to choose any normalized superposition of 
bound and continuum states as a target state. For example, choosing 

/ = iV 7 exp[-7(r-ro) 2 ]e ikor (29) 

as target state allows us to find a laser field that drives the particle to a predefined 
expectation value for position and momentum: (0y|f|0y) = r and ((f>f\p\(f)f) = k . 
This choice is of course not unique since, for example, any choice for the width 7 results 
in the same expectation values. 

2.5.2. Local operator Instead of using a non-local operator like the projection operator 
we can also employ a local operator for O, e.g., O = /(r). The most popular example 
for the local operator is the S function, 

J, = (¥(T)|*(r-r )|¥(T)> = |*(T,r )| 2 , (30) 

which maximizes (in the single-particle case) the density at point ro- The more density is 
squeezed into the point the higher the yield J\. The target by itself does not sound very 
physical but in practice it allows us to maximize the density distribution in a given region 
in space. In the next section we will show that a function different from the 5 function 
corresponds to a multi-objective optimization. In the practical implementation the 
5 function will be approximated by a narrow Gaussian. 

A different choice for driving the density towards a target density nf [similar to 
equation (j2"8l ] is 

min || VMT) - V%ll 2 = - {2 - 2 /* 0.(r,7>,(r)} , (31) 

where we assume normalization for n(T) and n/. The minimization corresponds to a 
maximization of 



J x = Jdv^n{v,T)n f {v). (32) 

2.5.3. Multi-objective target operators Within the formalism established, we are not 
restricted to a single objective. We can also employ a multi-objective target operator 
like 

d = J2PA- (33) 

3 

For example, the operators Oj can be projection operators for different excited states. 
The weights /3,- are chosen to balance the different objectives. If /5L- is chosen negative, 
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the optimization will try to minimize the expectation value of Oj, e.g., the occupation of 



a specific excited state. Combinations of projection and local operators are also possible. 
Note that the choice of a sum of projection operators has to be distinguished from the 
projection on the coherent superposition of these states. 

A local operator which is not a 5 function corresponds to a multi-objective 
optimization. This can be shown by considering the limit of infinitely many 5 target 
operators 



where the weight function (3(r) can be identified with an arbitrary local operator. 

2.5.4- Finite penalty versus complete controllability Introducing a (positive) penalty 
factor ol has the immediate consequence that a final target state occupation of 100% 
cannot be achieved. This can be proven within a few steps by achieving a contradiction: 
Assume that we have found the optimal field e opt (i) which drives the system from 
the initial state 4>i = ^(0) to the target state \&(T) =0/. According to equation ( l23i) the 
initial state for the Lagrange multiplier is then x{T) = <Pf- Since the two Hamiltonians 
of equations (f20l) and ( l22l) are the same, the time-evolution operators for and x{t) 
are identical: 



*(T) = U(T,t)U{t,0)*{0) = X {T) 
tf(t) = x {t). 

Inserting this finding in equation ( fl9l gives 

a k e k {t) = -lm(^(t)\fl k \^(t)) = 0, k = x } y } 



resulting in the statement that e opt (i) = or that 100% overlap cannot be achieved. 
Except for the trivial case e(t) = which presents a minimum for the functional if the 
initial and target state are orthogonal, (<f>i\^>f) = 0, we have to deal with the fact that 
100% occupation of the final state cannot be obtained with this kind of algorithm even 
if the system is completely controllable in principle. 

2.6. Algorithm to solve the control equations 

In the following we present two approaches to solve the standard optimal control problem 
[see equations (GO), (jHJ), and (jHJ)]. The first approach is an iterative solution [33] of the 
equations f|T9|) . ([20]) . fl22|) . and fl23|) . It provides the starting point for all extensions 
developed in this work. The second scheme applies to the case where the target operator 
is a projection operator. In this case a slightly faster algorithm can be deduced. 

As the word iterative already indicates, it will be necessary to solve the time- 
dependent Schrodinger equation more than once. Even with the present computational 




(34) 



(35) 
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resources this limits the application of the algorithms to relatively low-dimensional 
systems. 

2.6.1. Standard iterative scheme The control equations ( fT9i) . ( 1201) . ( 1221) . and ( 123]) can 
be solved as follows [33]. The scheme starts with propagating fa = \I/( )(0) forward 
in time. For the initial propagation we have to guess the laser field e^°^(t). For most 
of the cases we consider, the trivial initial guess e(°)(t) = is sufficient. However, 
sometimes the algorithm gets stuck in this solution. In this rule of thumb for 

the initial choice is that the forward propagation has to result in a large enough value 
for ||6^(T)|| 2 . A strong dc-field e (0) (t) = const often proves to be helpful. This part of 
the iteration can be expressed symbolically by 

stepO: *(°)(0) ^ *W(T). (36) 
After the initial propagation we determine the final state for the Lagrange multiplier 
wave function x {T) by applying the target operator to the final state of the wave 
function, 0\E^°)(T). The laser field for the backward propagation for x^ifyi ^ (t) is 
determined by 

f\t)= --lm( X ^m 3 \¥ k \t)), j = x,y. (37) 



a 



3 



The propagation from x^(T) to x^(T — dt) is done with the field e^(T), where we use 
X {T) and ty^(T) in equation (I3"7j) . The small error introduced here is compensated 
by choosing a sufficiently small time step. In parallel, we propagate ^^(T) backward 
with the previous field e^(t). This additional parallel propagation is only necessary 
if the storage of \l/( )(t) in the memory is not possible. For the next propagation step 
from X (0) {T-dt) to x m {T-2dt) we use ^(T-dt) and x {0) {T-dt) in equation (157} . 
We repeat these steps until x (0) is reached. To check the reliability of the parallel 
propagation we project \&( '(0) onto fa and compare with 1. We summarize the whole 
iteration step by 



stepk: [*( fc )(T) ^ *W(0)] 

0*(*)(T) = x (k) {T) x (fc) (0). 

The last part of the zeroth iteration step consists in setting ^(^(O) = fa and propagating 
^^(O) forward with the field e^\t) determined by 

ef +1 \t) = - -\m{x {k \m \¥ k+1) {t)), j = x,y, (39) 
«; 

which requires the input of x^(t). Again, we have to use the saved values from the 
backward propagation or propagate from x (0) to x(T) in parallel using the previously 
calculated field lf°\t). We end up having calculated e^(t) and ^f^(T) which can be 
expressed by 

[X (fc) (0) x W (T)} f4m 

& = (0) — T ^( fc+1 )(T). 
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This completes the zeroth iteration step. The loop is closed by continuing with equation 
fl38|) . i.e., propagating O^^fT) = x^^(T) with e^(t) [see equation (l3Tj) ] backwards to 
X «(0). 

If the initial guess for the laser field is appropriate the algorithm starts converging 
very rapidly and in a monotonic way, meaning that the value for the functional J in 
equation fflUj) is increasing at each iteration step. The monotonic convergence can be 
proven analytically [32] [33] - In the proof an infinitely accurate solution of the time- 
dependent Schrodinger equation is assumed. Since this is not possible in practice, it 
may happen that the functional decreases in the numerical scheme, e.g., when absorbing 
boundaries are employed. This sensitivity provides an additional check on the accuracy 
of the propagation. We can summarize the complete scheme by 

stcpO: *(°)(0) e< -^4' ) #(°)(T) 

stcpk: [*M(T) 6< ^4' ) *( fc )(0)] 

6#«(T) = X {k) (T) X {k) (0) 

[x «( ) ^ X «(T)] 

=*(*+!) (0) e< ^« (T). 

2.6.2. Projection operator - rapidly convergent scheme If the target operator is a 
projection operator O = \<f>f)(<f>f\, a scheme with even faster convergence can be derived 
from the modified functional 

J 3 [e,*,x] = -21m U^{T)\<f> f ) J dt (x(t) 

The scheme differs only in two points. First, the iteration is started with propagating 
x(t) backwards in time which is possible if we set |xCO) = \4>f)- Second, the equations 
which determine the laser field [equations ( 1371) and ( |39l) ] are replaced by 

f\t) = -Um{(¥ k \t)\ X ^(t)) ( X (k) m 3 \¥ k \t))}, 
a>j 

ef +1 \t) = - -Im{ (¥ k+1 \t)\ X W(t)) (x^m^^it)) } , j = x,y. 
We can summarize this scheme by 

stopO: / = x (°)(T) ^ X (0) (0) 

stopk: [ X «(0) X (k) (T)} 

[M,W( T ) #«(0)] 
0/=x ( fc+ i) (T) ^ <fc I^W x (*+i)( ). 

Again, monotonic convergence can be proven. The scheme contains some freedom in 
the choice of the first overlap (ty( k ^ (t)\ X ^ (t)} (with k' = k or k! = k + 1) in equations 
(USD and (@ID, since 

(*(*)lx(*)> = <*(Olx(0>- 

The authors of Ref. [32] report that the convergence changes for different choices of 
(ty (t)\ X (t)) . In our implementation we update this overlap at every point in time. 



id t -H(t)) . (42) 
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2.7. Example: Two-level system 

Let us apply the previously developed theory to a two-level system. A brief review 



about the theory of two-level systems can be found in Appendix A. 2 For the two-level 
system the integrals which determine the field, i.e., equations (|4"3I) and fpHl) . reduce to 

m)\x(t)) =g* a (t)h a (t)+g* b (t)h b (t), 

(x(Wm) ) = PabK^g.it)^-^ + PbaK^g^t)^-^ 

= /i (K(t)g b (t)e-^ + hlit^it)^) , (46) 

with Cfc(t) = (k\ty(t) ) and where guit) = Ck(t)e luJkt was defined in equations ( 1A.5I) and 
(1A.6[) . The coefficients of the Lagrange multiplier wave function x{t) in the basis \a) 
and |6) are (%(*) ) = l k (t) = h k (t)e~ iuJkt . 

Our goal is to find a laser pulse which transfers the ground state to the excited 
state | b) at T = 400. For this purpose, we use the algorithm described in equation ( )45l) 
with the penalty factor a = 1.0 and the initial guess e(t) = 0.05. After 5000 iterations 



we obtain an excited state occupation of 0.9996 and the laser field shown in figure l(a 



Note that the functional tries to find a laser field which produces a high occupation and 



has a low fluence. This behavior is clearly visible in figure 1(c), where the target yield 
Ji = \ (b\^/(T))\ 2 [(— — ) line] jumps to 0.9945 after the first iteration which corresponds 

to a fluence of E = 0.9204 [( ) line], while the converged field has the fluence 

of E = 0.0786. The optimal laser field [see figure l(a)| has a constant amplitude of 
A = 0.02 and frequency to = 0.1568 



J=\m(T))\ 2 -a fdt 6 2 (t), 
Jo 

where we have dropped the third term J3 since it is always zero. 

The optimal pulse within the validity of the rot at ing-wave- approximation (RWA) 
has to have a constant envelope A(t) = const. This can be understood with the help of 
the pulse-area theorem [50], 

li [ dt A(t) = Ti, (47) 



J 

which states that the inversion of a two-level system (within the RWA) is achieved if the 
area under the pulse envelope A(t) multiplied with /z, i.e., the dipole matrix element, 
becomes it. Now consider the following functional: 

L = J^dt A 2 (t) J^dt A(t) - tt \ . 

The variation of L with respect to A yields the pulse-area theorem [equation (T4T|) ]. while 
the variation with respect to the pulse shape A(t) results in 

2A(t) = Xfi. (48) 

If we plug this result back into equation fj47|) and solve for A we find 
2tt 



A = (49) 



Ait) - ± (50) 
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Figure 1. (Color online). We apply optimal control theory to invert the population 
of a two- level system, (a) Optimized laser field of the 5000th iteration with a = 1.0. 
(b) Time evolution of the occupation numbers for the system propagated with the 

optimized pulse. The ground state occupation corresponds to the ( ) line and 

the excited state to the ( ) line, (c) Convergence of J 1 [( ) line], J [( ) 

line], and the fluence E [( ) line], (d) J 1 (O ) and E Q (□) for different penalty 



factors after 2000 iterations. The filled symbols [J,: 
iterations. 



(• ), E Q : (■)] correspond to 100 



We may argue that the numerical algorithm has converged to the optimal pulse, because 
the RWA is perfectly valid for the chosen pulse length. 

We have run the optimization for different values a of the penalty factor. The 
occupation J\ (• ) and the fluence (■) after 100 and 2000 iterations (O , □) are shown 
in figure 1(d) We observe that 100 iterations are not enough to maximize the functional 
for small penalty factors. Although the occupation jumps to values J\ > 0.99 within a 
few iterations, the comparison with the longer iteration shows that there is still room for 
improvement. Iterating long enough yields a similar fluence for a range of a from 0.8 to 
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6.0. For penalties a > 2.0 the occupation starts to drop significantly because the fluence 
term is over-weighted. Selecting to small penalty factors leads to numerical instabilities 
which can be compensated by increasing the numerical accuracy of the propagation 
algorithm until the propagation becomes to costly. 

The applied optimization method does not provide a possibility to find optimal 
fields with a predefined fluence E . This can be achieved only indirectly by the penalty 
factor. However, the mapping between the penalty factor a need not to be invertible. A 
practicable way to find laser pulses with a given fluence is presented in section section I3TT1 

2.7.1. Short time transfer In the weak field regime, i.e., where the propagation time is 
usually long enough to justify the use of the RWA, the application of OCT for two-level 
systems does not seem to be appropriate due to the large numerical effort compared to 
the simple formula [equation f|47|) ]. 

In table Q] we show a comparison between the target state occupation achieved 
with pulses obtained from equation (14TI) and from OCT. The results show that for high 
requirements on the inversion efficiency (> 0.995) OCT becomes inevitable already for 
five-cycle pulses (T = 200). If one requires an inversion efficiency around 0.90 the RWA 
is appropriate up to single-cycle pulses (T = 40). 

The superiority of the OCT method for short pulses will become even more apparent 
for a system with more than two levels [51]. In this case a high strength of the field 
(oscillating with resonance frequency) will result in excitations to other levels which are 
minimized by the OCT pulses. 

Table 1. Comparison of the yield P = \(b\^(T))\ 2 when propagated with the laser 
field obtained with the two-level (RWA) estimate versus the laser field from OCT. Note 
that the period of the oscillation with ui ba is T p = 40.08. In all optimal control runs 
we set the number of iterations to 5000. In columns four and five we show the fluences 
calculated from the RWA and the optimized pulse, respectively. The penalty factor 
(in the last column) was chosen to give the best occupation for each optimization and 
a fluence comparable to the RWA values. 



T 


p 

RWA 


p 




UOCT 


penalty 


400 


0.9986 


0.9996 


0.0803 


0.0786 


1.0 


200 


0.9944 


0.9996 


0.1606 


0.1592 


0.5 


100 


0.9774 


0.9991 


0.3212 


0.3402 


0.3 


50 


0.9897 


0.9996 


0.6417 


0.7743 


0.3 


40 


0.8567 


0.9917 


0.8030 


0.7091 


0.3 


25 


0.7696 


0.9990 


1.2430 


1.4569 


0.3 



2.8. Example: Asymmetric double well 

In the remaining examples we will focus on a one-dimensional asymmetric double well 
to test our algorithms. The double well is similar to that in reference [52] but features 
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an additional cubic term: 

V(x) = ^-^ + (3x\ (51) 

with ujq corresponding to the classical frequency at the bottom of the well and the 
parameter B adjusting the barrier height. The number of pairs of states below the 
barrier is approximately B/ujq. Here, we choose B — Uq — 1.0 and (3 = 1/256 which 
leads to two states below the barrier, as shown in figure El In order to analyze the laser 



0.5 



7 o 

>. 

S -0.5 



-1 

-6 -4 -2 2 4 6 
x [a.u.] 

Figure 2. The plot shows the model potential with the ground state ( ), the 

first excited state ( ), the second excited state ( ) and the third excited state 

( — • — ). Each state is shifted according to its eigenvalue. 

pulses from the optimization runs we calculate the excitation energies (see table [2]) and 
dipole moments (see table [3]) of the system by propagating in imaginary time. 

Table 2. Excitation energies in atomic units [a.u.] for the ID asymmetric double well, 
calculated by imaginary time propagation. 

10) |1) |2) |3) 
10} 0. 

|1) 0.1568 0. 

|2) 0.7022 0.5454 0. 

|3> 1.0147 0.8580 0.3125 0. 

|4) 1.5294 1.3726 0.8273 0.5147 




The time-dependent Schrodinger equation for the ID double well is solved on an 
equidistant grid, where the infinitesimal time-evolution operator is approximated by the 
2nd-order split-operator (SPO) technique [53] : 

/ rt+At \ 
fjt+At = Texp / _ • ^ dt , 

w exp(~f At)exp{-iV{t) At)exp(-^f At) + 0{At 3 ). 

Following the rapidly convergent scheme described in section 12. 6[ one needs four 
propagations per iteration (if we want to avoid storing the wave function). Within 
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Table 3. Dipole matrix elements for the ID asymmetric double well, calculated by 
imaginary time propagation. 





10) 


|1> 


|2) 


|3) 


|4) 


10) 


-2.5676 










11) 


0.3921 


2.3242 








|2> 


0.6382 


-0.7037 


-0.5988 






|3) 


-0.3865 


-0.4630 


1.7051 


0.1958 




|4> 


-0.1414 


0.2118 


0.1593 


-1.7862 


-0.0939 



the 2nd order split-operator scheme each time step requires 4 Fast Fourier Transforms 
(FFT) [51] for the backward propagations, because we have to know the wave-function 
and the Lagrange multiplier in real space at every time-step to be able to evaluate the 
field from equations ( ]4"3l) and (jHj). This sums up to 16 FFTs per time step and iteration. 



Table 4. Employed numerical parameters (atomic units). 



parameter 


T 


400.0 


pulse length 


X 

max 


30 


grid size 


dx 


0.1172 


grid spacing 


dt 


0.001 


time step 


e (0) 


-0.2 


initial guess 




10~ 5 


convergence threshold 



The parameters used in the runs are summarized in table HI The initial guess for 
the laser field was e^°'(t) = —0.2. This choice is arbitrary but has the advantage of 
producing a significant occupation in the target state at the end of the pulse, necessary 
to get the iteration working. Although the simple choice e°)(i) = 0.0 will work as well 
in most cases, it represents a minimum of the functional since initial and target state 
are orthonormal. Therefore the algorithm could get stuck in principle. 

In the following we apply the rapidly convergent algorithm to find the optimal field 
that transfers the ground-state to the lst-excited state. Choosing the penalty factor 
a = 2.2 the algorithm converges after 515 iterations to J = 0.8470. We consider the 
value as converged if the change of the functional between two subsequent iterations is 
smaller than AJ fe_1 ' fc = 10 -5 . 



The optimal laser field is shown in the upper panel of figure 3(a) Applying this 



laser field to the system yields an occupation of 0.9944 in the target state. The laser 
field exhibits a fluence of 0.0670 which is 16% less than a monochromatic pulse with a 
similar final occupation would need. The first step to analyze the optimal pulse is via its 
spectrum shown in figure 3(b) It is dominated by three narrow peaks which are located 
close to the field-free excitation energies u i = 0.1568, = 0.5454, and U02 = 0.7022. 
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Further we observe a group of peaks around u> i (u e [0,0.1] and u G [0.2,0.3]) which 
do not correspond to any excitation energy of the field-free Hamiltonian. However, 
these frequencies play an important role in the transition process. If we filter out these 
frequency components, rescale the fluence to E = 0.0670, and then propagate this 
modified laser pulse, we find at the end of the pulse the following occupations: ground- 
state 17.0%, first excited state 24.4%, second excited state 58.4%, and in all higher levels 
0.2%. In particular, the direct transition and the back transfer from the intermediate 
level 1 2) to the target state in the indirect process are less efficient without these extra 
frequencies. Further analysis of this kind shows that the low-frequency components 
and especially the zero frequency component (bias) are crucial since they introduce a 
(slight) shift compared to the (field-free) resonance frequencies, visible as a broadening 
of the uqi peak in figure 3(b) On the other hand, if these components are missing the 
remaining frequencies become (in that sense) off- resonant, resulting in a low occupation 
of 24.0% of the first excited state. 

If we filter out everything except for the extra peaks we find a target state 
occupation of 5%. Understanding these extra peaks as a third type of transfer process 
suggests that a mixing of transition processes in this case seems to be superior in terms 
of the maximum target yield per fluence than a simple monochromatic pulse. 

The gain in the occupation of 0.01% and in the fluence of 16.25% compared to 
a simple monochromatic pulse has a high price - the optimized pulse is much more 
difficult to realize in an experiment. Although the gain improves with shorter pulse 
lengths (see tabled]), this example demonstrates the typical dilemma between theory 
and experiment: Calculated pulses often have a far too complicated spectrum to be 
produced in practice. In sections 13.41 - 13.51 we demonstrate how this dilemma can be 
solved. 

To conclude the analysis we look at the convergence behavior of the applied scheme 
[see inset of figure 3(b)| . We find a fast convergence within the first 20 iterations. After 
these iterations the improvement of the yield slows down. 



3. Constraints on the optimal laser field 

Despite its importance only a few attempts have been made to take further restrictions 
on the optimal field into account. In Ref. [29] a scheme to calculate the pulse for a given 
fluence is shown. However, it does not make use of the immediate feedback introduced 
in Ref. [32], and it suffers from a rather unstable convergence. A constraint on the 
spectrum is considered in Ref. [55] for a steepest descent method which, in the quantum 
control context, also exhibits a poor convergence and a strong dependence on the initial 
pulse [56]. An elegant way to restrict the spectrum has been presented in Ref. [57] . 
This scheme preserves the rapid and monotonic convergence behavior of the underlying 
scheme [32] by projecting out those parts of the time- dependent wave function which 
are responsible for the unwanted spectral components. However, this method is not 
sufficiently general and does not easily allow for an additional fluence constraint (as it 



Quantum Optimal Control Theory 



18 




Figure 3. (Color online). Optimization of the |0) — » |1) transition, (a) Top: 
Optimized field. Bottom: Time evolution of the occupation numbers | (i) | xz.) | 2 [n = 

( ), n = 1 ( ), n = 2 ( ), and n = 3 ( — • — )]. (b) Spectrum of optimized 

field. The vertical lines indicate the transition frequencies lo 01 , w 12 , and lo Q2 - Inset: 

Convergence of J l [( ) line], functional J [( ) line], and the fluence E Q [( ) 

line; scale on the right]. 



keeps the penalty factor). 

The schemes shown in the following are similar to what we have discussed in 
Ref. [51], but are presented in a slightly more general way This allows us to incorporate a 
large variety of experimental constraints in the optimization, for example, fluence and/or 
spectral constraints and phase-only shaping. As we will see in section [37^1 and section |3~5| 
they show very good convergence, although a proof of monotonic convergence similar to 
reference [33] is not possible here. The difficulty consists in atj which is changing during 
the iteration and in the case of spectral constraints due to the (brute-force) modification 
of the field. But, even for the "brute-force" spectral filter we find a good convergence 
unless not too many essential features of the pulse are suppressed. 

Since we do not expect a monotonic convergence we have to add some additional 
intelligence to the algorithm, i.e., we store the field which produces the pulse with the 
highest yield in the memory. This field is considered as the result of the optimization. 



3.1. Fluence constraint 



In order to fix the fluence of the optimized laser pulse to a given value E , we have to 
replace the functional J2 by 

f -T 



E 



alt ej(t) - E 



(52) 



Here otj is a (time-independent) Lagrange multiplier. Instead of specifying atj, we have 
to prescribe specific values E 0j for the components E 0x and E 0y of the laser fluence. 
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Since (here) atj is a Lagrange multiplier we have to vary with respect to it when 
calculating the total variation of J [cf. equation fTTTj) ]. The variation with respect to a k 
results in an additional equation 



alt e k \i) = E { 



k 



x,y. 



(53) 



In the case where a k is a penalty factor its value has to be set externally, while here 
the additional equation can be rewritten [29] to determine the value of the Lagrange 
multiplier a k . Inserting equation ffl9l into equation fl53l yields 



or 



Oi k 



T 



dt 



lm( X (t)\jl k \*(t)) 



~-w k {t) 



E, 



k 



Ok 



(54) 



where fi k is the dipole moment operator. 

The remaining part of the functional stays the same, so the variations do not change, 
and we keep the control equations: (fT9l) . (120]) . ff22l . and ff23l . 



3.1.1. Algorithm The set of coupled equations which have to be solved is now given 
by the equations ( fT9l) . ( 120|) . ( |22|) . ( 1231) . and ( 1541) . The scheme below shows the order in 
which these equations are solved in the /cth step. 

e( fe )(t) 



step k: *( fc )(0) 



[*(*)(T) 



X (*)(T) 



e< fc )(t) 
e«(t) 



*W(0)] 

x (fc) (o), 



with the laser fields e^(t), e^ k \t) given by 



(A:) 



Im<X«(t)|fc|*«(t)>, 



,(fc+i) 



a 



(fc+i) 



J — x iVi 



where the Lagrange multiplier aq fe is defined by 



(fc+i) 
or- 



\ 



Io dt 



(k)~(k) 



-i 2 



J =x,y. 



The initial conditions in every iteration step are 
tf(r,0) = &(r), 
X (r,T) = 0tt(r,T). 



(55) 



(56) 
(57) 



(5f 



(59) 



The scheme starts with the propagation of *ff(°'(t) forward in time using the laser field 
e^(t) which has to be guessed. The result of the propagation is the wave function 



Quantum Optimal Control Theory 



20 



\I>( )(T) which is now used to calculate x {T) by applying the target operator [equation 
([59]) ]. We continue by propagating x^(t) backwards in time using the laser field e^(t) 
defined by equation ( 1561) . To solve equation ( 1561) . we have to know both wave functions 
tf(°>(t) and X {0) (t) at the same time t, which makes it necessary to either store the 
whole time-dependent wave function \l/( )(t) or propagate it backwards with the previous 
laser field The avoided storage is indicated by the brackets in the scheme ( 155]) . 

Moreover, it is necessary to provide an initial value for at,- which we choose to be 



(o) 
a- 

3 



& <f (*) 



(o). 



jo [ j 



3 =x,y. 



The result of the backward propagation x^ (t) is the laser field e^ ^ (t) which is rescaled 
to the right value with equation ( 1571) giving e^(£). This completes the first iteration 
step. The second (k = 1) or, in general, the kth step repeat the described procedure 
starting with the initial state ^^(O) = fa and the rescaled field e^ k \t). 

The scheme described above has some aspects in common with the techniques 
described in Refs. [29], [32]. The basic idea of incorporating fluence constraints in the 
optimization algorithm was given in Ref . [29] . However, in contrast to this reference we 
make use of immediate feedback [see equation (|56j) ]. i.e., the backward propagation is 
accomplished by updating x(t) an d e(£) in a self-consistent way, which was suggested 
in Ref. [32] . On the other hand, the technique presented by the authors of Ref. [32] 
does not allow to build in fluence constraints, since aij is not a Lagrange multiplier in 
their case. Roughly speaking, the technique presented above is a combination of these 
approaches. 



3.2. Generalized filtering technique 

In contrast to the fluence constraint, the following general technique is not derived from 
a functional. Rather, a general filter is applied "brute force" in every iteration to the 
laser field, 

eoutjO) = fen j Wi- 
lli principle, the filter Q can be any operator. A few examples are discussed in sections 

Since the functional itself stays the same as the standard functional [equations (J7J)- 
(TTOj) ]. we have to solve the usual set of control equations: (1191) . (1201) . ( 1221) . and (BHD - 



3.2.1. Algorithm The algorithm with built in general filtering is similar to the one 
presented in the previous section (section 13.11) with two important differences: The 
factor ay is a penalty factor. It has to be specified from the start and remains unchanged 
during the optimization. Second, the update of e^ 1 ^ (t) in equation (1571) is replaced by 



Sk+i) 



(t)=G ~ef\t) 



J =x,y, 



(60) 



where the symbol Q indicates a given filter operator acting on the field €j(t). 
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3.2.2. Spectral constraints If spectral filtering is required, we formulate the constraint 
with the help of a filter function fj(oo), the Fourier transform J 7 , and its inverse 
Equation ( 1601) is now replaced by 



ef +1 \t) = T' 1 [fM? {t)\\ , J = x,V. (61) 

Since 6j(t) is real- valued we have to make sure that fj(w) = fj(—u). For example, the 
filter function could be chosen to be 

fj(u) = exp[-7(^ - uj ) 2 } + exp[-7(u; + ^ ) 2 ], (62) 

so that only the components around the center frequency ±ujq are allowed in the pulse. 
If one uses instead 

fj(u) = 1 - (exp[-7(u; - cu ) 2 ] + exp[->y(uj + u; ) 2 ]) , (63) 

one would allow every spectral component in the laser field except the components 
around ±u> . 

In Ref. |58j we show that the spectral filter technique can be derived from a modified 
form of the standard functional. 



3.2.3. Laser- envelope constraints Even though the formulation of the theory with a 
time-dependent penalty factor provides already one way to enforce a time-dependent 
shape function hj(t), we want to present an alternative way. Here, we replace equation 
(EOD by 

ef +1 \t) = h J (t)~ef\t), j = x,y. (64) 

While the more elegant way is to use a time- dependent penalty factor [48J , its application 
is not always possible, for example, in the case of a fixed laser fluence where the penalty 
factor does not exist. While the alternative method still enables us to impose restrictions 
on the laser envelope and at the same time to fix the fluence to a given value. 



3.2.4- Phase-only shaping Many experiments are carried out using only phase-shaping, 
only the phases of the spectral components are optimized but the amplitude 



i.e. 



spectrum itself stays fixed (for details on the experiment see section Appendix A.l). 
This is done to reduce the enormous search space for the genetic algorithm and to 
achieve a faster convergence. 

For the implementation of phase-only shaping into the computational optimization 
we have to replace equation (I60l) by the following equations: 



3 



(t) 



,(*+!) 



(t) = T- 



(*)/ 



3 = x,y, 



(65) 
(66) 



, ^)l 

contains the predefined amplitude spectrum which, in the 
experiment, corresponds to the spectrum of the laser pulse that enters the pulse shaping 
device. 



where the function Aj(u 
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3.2.5. Combination of filters In principle, the filters can be freely combined. Then 
equation (1601 has to be extended to 



,(*+!) 



(*) = Q. 



N 



<N-1 



...Qx 



3 



(*) 



J = x,y. 



(67) 



Care must be taken if the filters are conjugated, for example, a filter in the frequency 
domain will change the field in the time domain as well. In this case the order of the 
filters is important. For example, we want to restrict the spectrum and at the same 
time require a Gaussian shaped laser envelope. In general, it is impossible to satisfy 
both requirements at the same time. However, if the spectral filter function is broad 
enough (unlike a 5 function which will always result in a constant envelope) the time 
filter can be satisfied in a reasonable way. Of course, if the filters are conjugated only 
the last filter will have the full desired effect. 



3.3. Generalized filtering technique with fluence constraint 

It is also possible to combine the fluence constraints with the filtering techniques. This 
combination makes it possible to implement even more realistic experimental constraints 
in computational pulse optimizations. Therefore we use the functional and the control 
equations discussed in section 13.11 and add the "brute force" method of section 13.21 



3.3.1. Algorithm The algorithm is a simple combination of the schemes presented 
above. Basically we use the scheme (1551) . equation (1551) . and instead of equation ( 1571) we 
employ 



T(t) 



,(fc+i) 



(*) 



a 



'3 
(k) 



' -A *(*)■ 



a 



(fc+i) 



J = x,y, 



(68) 
(69) 



where Q is a given filter operator and otf +l ^ is evaluated inserting the filtered field ej fe ' (t) 



m 



(fc+i) 
Or- 
3 



\ 



ftdt af¥«\t) 



Eo 4 



J = x,y, 



(70) 



to enforce the predefined value for Eo r Note that the total spectral power is related to 
the time-integrated quantity by Parseval's theorem, 

dt 6(t)6(T - t) [e^t)] 2 = — dcu |e»| 2 . 

-oo ^ J ~ oo 

In this combined form of the algorithm we first apply the filter function to the laser field 
( 168]) and then rescale the field to yield the right value for E . [see equation (169]) ]. 
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3.4- Example: Direct Transition 



We now present the results of the algorithm with spectral constraints and a penalty 
factor. In the previous chapter the transfer of the particle occured via a mixture of a 
direct transition and indirect transitions. This motivates the following aim: Find a laser 
pulse that produces a high yield and that contains only spectral components centered 
around the resonance frequency cjqi- 

In order to find an optimal pulse with ujqi as the center frequency, we use a Gaussian 
frequency filter f(u) according to equation (|62|) centered around u = u i and with 
7 = 500. 

After 50 iterations the algorithm finds a laser pulse which yields of 99.97%. The 
penalty factor has been set to a = 0.05. The obtained value for E = 0.090 is slightly 
higher than the estimate from the two-level model (Eq = 0.08, J\ = 99.30%). The 



slight envelope on the field, visible in the upper panel of figure 4(a), is caused by the 

finite width of the Gaussian [see ( ) line in figure 4(b)] . Frequency components 

near c^oi are still allowed in the pulse and result in a beat pattern (envelope). The 
time-dependent occupation numbers confirm that the higher states are not occupied 
during the transition [see the lower panel of figure 4(b)] . The convergence shown in 



figure 4(b) is rather smooth. Note that if we desire a sinusoidal field with a constant 
envelope, we have to reduce the width of the Gaussian to a 5 function to allow only one 
single component in the spectrum. Using such a filter we obtain a yield of 99.79% and 
Eq = 0.085. This field oscillates with the amplitude A = 0.0207 which is slightly higher 
than the amplitude derived from the pulse area theorem [see equation ( lA.li 
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Figure 4. (Color online). Optimization of the |0) — > |1) transition with a Gaussian 
frequency filter around u Q1 — 0.1568. (a) Upper panel: Optimized field. Lower panel: 

The time-dependent occupation numbers confirm that only the ground state ( ) and 

the first excited state ( ) take part in the transition process. The second excited 

state population ( ) is hardly visible, (c) Spectrum ( ) and filter function f(uj) 

( ), scaled by O.Of. Inset: Convergence of J 1 . The (O ) indicates the iteration 

with the highest yield. 
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Figure 5. (Color online). We apply the optimization algorithm for the transition 
|0 ) — > |1 ) with a double Gaussian frequency filter allowing only frequencies around tu Q3 
and oj 31 , and in addition we set E Q = 0.320. (a) Optimized field and time-dependent 

occupation numbers [ground state ( ), first excited state ( ), and second excited 

state ( )]. (b) Spectrum and filter function f(uj) ( ), scaled by O.Of. Inset: 

Convergence of J v The (O ) indicates the iteration with the highest yield. 



The process |0) — > |3) — > |1) using the third excited state as an intermediate state 
plays only a minor role in the examples considered above. Now, we are going to optimize 
the laser pulse such that the transition occurs exclusively via this process. In addition 
we require E = 0.320. This time a double Gaussian filter is centered at co> 13 and u 03 . 
The width parameter is again 7 = 500. 

The results are shown in figure [5j Like in the previous example, the high restrictions 
within the optimization lead to a rather erratic convergence [see inset of figure 5(b)| . 
The field, shown in the upper panel of figure 5(a) , corresponds to the 162th iteration and 
produces a target state occupation of 99.89%. The time-dependent occupation numbers 
[see lower panel of figure 5(a)| show that the transition occurs exactly in the desired 
way. 



4. Time-dependent control targets 

The control targets we have considered so far refer to the maximization of a quantity 
at the end of the laser pulse, e.g., the occupation of some excited state. This kind of 
optimization objective is called a time-independent target, since it leaves the dynamical 
path the quantum system follows towards a target state undefined. In this chapter we 
demonstrate that it is also possible to control this path, i.e., to find the laser pulse which 
leads the quantum system as close as possible along a predefined trajectory [59]. The 
path could be simply a trajectory in the configuration space but it may also be a path 
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in a more abstract sense, e.g., a trajectory in quantum number space to control how a 
transition takes place. 

Control targets that require a time- dependent formalism are the control of bond 
distances in molecules (e.g., steering the fragmentation process in time or using a laser 
to keep a certain bond distance), the optimization of high harmonics [HI EH], and the 
control of currents in time, e.g., in a molecular switch. 

To our knowledge, three different methods for the control of time-dependent targets 
have been proposed so far: A 4th-order Euler-Lagrange equation to determine the 
envelope of the control- field has been derived in Ref. [60]. However, it is restricted 
to very simple quantum systems. 

A very elegant method, known as tracking has been proposed in Refs. [61] [62] . 
Despite its tremendous success, this method bears an intrinsic difficulty: One has to 
prescribe a path that is actually achievable with a laser field, otherwise singularities in 
the field appear due to the one-to-one correspondence between the laser field and the 
given trajectory. In practice, this may require a lot of intuition. However, the most 
severe drawback is that tracking cannot be used together with constraints on the laser 
pulse (see section [3]) , which is again a problem of the one-to-one correspondence. 

The third method is an optimal control scheme for time-dependent targets [63] 164"]. 
The new method is monotonically convergent and in contrast to tracking it does not 
rely on one's intuition in choosing the right control targets. Furthermore, the method 
is not restricted to two-level systems and can be extended to incorporate fluence and 
spectral constraints, as we will show in section 14.51 First applications of the method 
can be found in [631 ESJ M> ESI M\ ■ 

In section H] we present the formalism for the control of time-dependent targets [59] 
in the same way as in section [2721 The algorithm [63], [66] to solve the resulting equations 
is described in section 14.41 

In section 14.51 we present a novel algorithm which combines the standard control 
algorithm of section 14.41 with the constraints discussed in section [3j We test the 
developed algorithm for the asymmetric double well model (see section 12.81) and discuss 
the results for the control of a path in quantum number space in section 14.61 



4-1. Derivation of the control equations 

Consider a modified form for the first term of the standard functional in equation (1101) 

jm = ^ j o T dt w(t)(*(t)\6(t)m)), (71) 

with w(t) representing a time- dependent weight function which is normalized in the 
following way 

i-T 



[ dt w(t) = T, w(t) > V t. 
Jo 



The weight function w(t) is supposed to steer the relative importance of the time- 
dependent target operator. If the target operator 0{t) is positive semidefinite then J\ 
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will reach its maximum, if at each point in time the expectation value of the operator 
(fy(t)\0(t)\ty(t)) is maximized. For the moment we will keep the operator as general as 
possible and postpone the discussion of different examples to section 14.31 

The other parts of the standard functional, J2 in equation (jSj) and J3 in equation 
([9]), remain unchanged. Thus, the functional derivative of J with respect to \l/ becomes 

<)Jl kw(r)d(r)^*(r',T), 



which yields 



<5^(r',r) T 

^T^y = - i fa + H(r)) X*(r', r) - [**(„, t)5(t - r)] ^ 



S„J = £dr ^w(t) (*{t)\6\6*(t)) + i (fa - H(rj) x(t)|**(t))} 
-(X(31W) + (X(0)|«*(0)). (72) 



=0 



4-2. Time- dependent control equations 

Setting the variation with respect to ^ in equation (1721) equal to zero, we obtain an 
inhomogeneous time-dependent Schrodinger equation for the Lagrange multiplier x(r, t) 

fa - H{t)) X (r, t)= - l W (t)d(t)*(r, t) , X (r, T) = 0. (73) 
Its solution can be formally written as 

X (r, t) = Ulx(r, to) - ^ jf dr # far)6(r) *(p, r)) , (74) 
where ?7* is the time-evolution operator [68] defined by 



E/J=Texp 



dt' H(t') 



to 



Since the other parts of the functional correspond to the standard functional, also the 
variations with respect to x(r, t) and e k (t) are identical to equations (TSUj) and (JEJ), which 
we restate here for convenience 

fa - £(t)) *(r, t) = 0, *(r, 0) = &(r), (75) 

a k e k (t) = - Im(x(t)|Afc|*(*)>, fc = a; > 3/. (76) 

Similar to section 13.11 we can introduce fluence constraints by using J 2 [equation (152]) ] 
instead of J2 [equation (jSJ)], for which we obtain an additional equation from the variation 
with respect to the Lagrange multiplier a k 

[ dt e k 2 {t) = E 0k . (77) 
Jo 

The set of equations that we need to solve is now complete: (173|) . (1751) . (1761) . and in the 
case of a predefined fluence, we have in addition equation (1771) . 
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4-3. Target operators 

The physical meaning of the functional J± given by equation ( J7T1) depends on the choice 
of the target operator. In the following we present the most important choices and 
discuss their physical interpretation. 



Final-time control: Since our approach is a generalization of the standard optimal 
control formulation given in section 12.21 we first observe that the latter is trivially 
recovered as a limiting case by setting 

w(t) = 2T5(t-T), 6 = \<f> f )((f) f \, 

where we use the definition: J Q dt 5(t — T) = 1/2. Here \<f>f ) represents the target state, 
which the propagated wave function ^/(t) is supposed to reach at time T. In this case 
the target functional reduces to [29l [32] 

</i = i<*cni<Mi 2 . 

The target operator may also be local, as pointed out in Ref. [33] . If we choose 
w(t) = 2T5(t — T) and O = 5(r — ro) (the density operator), we can maximize the 
probability density in r at t = T, 

Ji = <*(T)|6|*(T)) = n(r ,T). (78) 

In the actual calculations, the 5 function can be approximated by a narrow Gaussian. 



Maximizing the average: In the literature [60] 163] the functional (ITT]) has so far only 
been used with a time-independent target operator, e.g., 

o = |0/>(0/|, 

combined either with a time- independent (w(t) = 1) or time-dependent weight function 
[69]. 

In the first case, we require the laser pulse to maximize the average occupation in 
state i.e., the earlier the laser pulse drives the initial to the target state (and keeps 
it there), the higher the yield. 

Wave function follower: The most ambitious goal is to find the pulse that forces the 
system to follow a predefined wave function <p(r,t). If we choose 

w(t) = 1, 

6{t) = m))m)\, (79) 

the maximization of the time- averaged expectation value of 0(t) becomes almost 
equivalent to the inversion of the time-dependent Schrodinger equation, i.e., for a given 
function 0(r, t) we find the field e(t) so that the propagated wave function \l/(r, t) comes 
as close as possible to the target 0(r, t) in the space of admissible control fields. We can 
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apply this method to the control of time-dependent occupation numbers, if we choose 
the time-dependent target to be 

|0(t)) = c (t)e-^|0> + Cl (t)e- i£l *|l) + c 2 (t)e- i£2t \2) + . . . , (80) 
H \n) = £ n \n) 

6{t) =W))W)l (si) 

The functions |c (t)| 2 , |ci(t)| 2 , |c 2 (t)| 2 , . . . are the predefined time-dependent level- 
occupations which the optimal laser pulse will try to achieve. In general the functions 
c (i), Ci(t), c 2 (t), . . . can be complex, but as demonstrated in Ref. [6_Ij, real functions can 
be sufficient to control the occupations in time. For example, if in a two-level system the 
occupation is supposed to oscillate with frequency Q we could choose co(t) = cos(fii) 
and, by normalization, Ci(t) = sin(f2t). This defines a time-dependent target operator 
by equations ( IHUi) and ( |8B . 



Moving density: The operator used in equation (1781) can be generalized to 

6(t)=«J(r-r (t)), (82) 
1 ' T 



Ji = f I dtw(t)(V(t)\5(r-T (t))\*(t)) 



T 



= - jl dt w(t)n(r Q (t),t). (83) 

Here J\ is maximized if the field focusses the density at each point in time in r (t). An 
example can be found in Ref. [67J. The relative importance of different time intervals 
can then be adjusted by the weight function w(t). 

4-4- Algorithm for time-dependent control targets 

Equipped with the control equations (1731 . (1751) . (175]) we have to find an algorithm to 
solve these equations for e(t). In the following we describe such a scheme which is similar 
to those in Refs. [311 E3] where the additional parameters rj and £ have been "artificially" 
introduced (not derived by a functional variation) to "fine tune" the convergence of the 
algorithm. A monotonic convergence in J can be proven if r] G [0, 1] and £ G [0, 2] [34]. 
Here, ex. is always a penalty factor. 

The algorithm starts with propagating ^( )(0) = <f>i forward in time with an initial 
guess for the laser field e(i)^\ 

stepO: *(°)(0) ^ *W(T). (84) 

The backward propagation of X (f) ls started from X^(T) = solving an 
inhomogeneous time-dependent Schrodinger equation which requires ^^(t) as input 
(with k = 0), 

stepk: [*( fc )(T) ^ ^«(0)] 
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The brackets indicate that the storage of the wave function \&(°)(t) can be avoided if 
we propagate it backwards in time as well using e^(t). The backward propagation of 
X (t) requires the laser field determined by (k = 0), 

f\t) = (1- V )ef\t) -^-lm(x (k \m,\¥ k \t)), j = x,y. (86) 



The next step is to start a forward propag ation of *W(0) = & (k = 0) 

k«(o) ^ fc)(t) x « ( r)] 

[*(*)( ) ^ *W(T)] (8T) 

(0) £< —^ t] #( fc+1 )(T). 

and calculate the laser field 

ef +1 \t) = (1-C)i k \t) -^lm(x (k \tm^ k+1 \t)), j = x,y. (88) 

a>j 

If we want to avoid storing x^(£) i n the memory we have to perform an additional 
forward propagation which in turn requires the knowledge of \l/( )(t). These extra 
propagations are indicated by the expressions in brackets. After the time evolution 
is complete we can close the loop and continue with equation fl85l) . 
The whole scheme can be depicted symbolically by: 

stepO: *W(0) £< -^.' ) *(°>(T) 



stopk: [^( fc '(T) £< ^* ) *( fc )(0)] 

X (,) (T) ^C^C*) x(fc)(Q) 



[X (,) (0 ) ^«_^>M xW(T)] 
[tf«(0) e< ^' *W(T)] 
*(fe+i)(0) > w *( fe+1 )(T). 

The main difference between this iteration and the schemes used in Ref. [31] is that 
one needs to know the time-dependent wave function ^(t) to solve the inhomogeneous 
equation (1731) for the Lagrange multiplier x{t)- 

The choice of rj and £ completes the algorithm. £ = 1 and r] — 1 correspond to 
the algorithm suggested in Ref. [33], while the choice £ = 1 and 77 = is analogous 
to the method used in Ref. [29] with a direct feedback of iff( k '(t). Further choices are 
discussed in Ref. [31]. A more detailed discussion on the convergence of this algorithm 
with exclusively time-dependent target operators and a modified version of the target 
functional can be found in Ref. |66|. 



(89) 



4-5. Algorithms for time- dependent control targets with constraints 

For time-dependent targets a spectral restriction of the laser pulse turns out to be even 
more important than in the time-independent case. Besides the need to incorporate 
experimental limitations, spectral restrictions become important already at the level of 
modelling. For instance, assume that we want to optimize a time-dependent occupation 
or density such that it oscillates at a given frequency but the optimized pulse is not 
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allowed to contain this frequency, e.g., to optimize high-harmonic generation, in that 
case a spectral constraint becomes inevitable. 

In the following we present an extension of the algorithm discussed in section [3] 
to incorporate time- dependent targets. During the presentation in section [3] we have 
distinguished between targets with a constraint on the ffuence, filter algorithms using a 
penalty factor, and fluence constraints together with filter algorithms. Here we combine 
all of them in one scheme: 



step k: ^(0) 



eW(() 



#(*)(T) 
[*(*)(T) 



X {k \T) 



e<»(t) 
eW(t),tf( fe )(t) 



tf(*)(0)] 
X «(0), 



(90) 



with the laser fields e^ k \t), e^ k \t) given by 

1 



fa) 



(i) 



(fe) 

i 



Im< x «(t)|,a,-|*M(i)}, 



3"W 



,(fc+i) 



(fe) 



(fe+1) 

j 



3 — x iVi 



(91) 
(92) 
(93) 



where ^ is a generic operator. Its precise form depends on the application, several 
examples can be found in section 13.21 In the case of a fixed fluence the Lagrange 
multiplier a,- fe is defined by 



a 



(fc+i) 



Io dt 



(fe)-(fe) 



(94) 



a,- 



otherwise (when is a penalty factor) we set af +1 ^ 
The initial conditions in every iteration step are 

tt(r,0) =&(r), 
x(r,T) = 0. 

If only the fluence has to be kept fixed to a given value we use £/[e(t)] = e(t). 



(95) 



4-6. Example: Optimal Control of time- dependent occupation numbers 

In this section we present an application of the algorithms fl89l) and fl90|) . We want 
to control the occupation numbers of the asymmetric double well model in time. We 
look for a laser pulse which drives the system as close as possible along a given path 
in quantum number space, defined as follows: First, excite the system from the ground 
state to the fourth excited state, then dump the occupation as fast as possible to the 
third excited state [realized by the step function 0(5 T/8 — t) in equation ( llOOp ]. and 
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finally transfer the occupation to the first excited state. In mathematical terms, the 
target occupation numbers |c n (t)| 2 are defined as follows: 

c<j(t) = 0(T/2 - t) cos(7rf/T) (96) 

Cl (t) = 6{t - 3T/4) sin(27rt/T - 3vr/2) (97) 

c 2 (t) = (98) 

c 3 (t) = 6(t - T/2) (1 - |c (t)| 2 - | Cl (i)| 2 - |c 4 (t)| 2 ) 1/2 (99) 

c 4 (t) = 0(T/2 - t) sin(vrt/T) + 9(t - T/2) 6(5 T/8 - t). (100) 

For t < the system is in the ground state. The choice of the timings and the 
particular functional form of the population curves are just examples to demonstrate 
the capabilities of the optimization method. 

The squared moduli of the functions in equations (196]) - (11001) are represented in the 



middle panel of figure 6(c) by the ( ) lines. The coefficients c n (t) define the target 

operator 

6 =W))W)\, 

4 

W)) = E c «(*) e_i£rif i n >- 

n=0 

The control objective has some resemblances with the one shown in section 13.51 where 
we have used a filter function to achieve a similar path-selective excitation. The time- 
dependent target operator, however, exhibits a much more precise control. We can 
specify not only the timings of the single transitions but also the rate. On the other hand, 
the choice of the trajectory does not require the same amount of intuition compared 
to the tracking method [61], i.e. with the help of OCT we can find a laser field that 
guides the system as close as possible along the trajectory, but not necessarily exactly. 
Therefore we are able to use also "unphysical" trajectories [M] , e.g., the step-like change 
at t = 5 T/8 in equation (11001) . where tracking methods will not work at all. 

In the following we present the results for the time-dependent target described 
above (see equations (1961) - (HOOD ). To suppress unwanted frequency contributions and 
to obtain a simple pulse shape we apply the filter algorithm described in section 14.51 
with the parameters a = 0.2, w(t) = 1, T = 800, e (0) (t) = 0, (7 = 1, rj = 1). The filter 
functions are Gaussian-shaped and centered around the field-free transition frequencies 
co> 04 = 1.5294, W43 = 0.5147, and U31 = 0.8580. In addition, we want to improve the final 
occupation of the first excited state and to decrease the importance of the fast transition 



at t = 500 by using a weight function which is shown in the top panel of figure 6(c) and 
defined by 



f(t) = 1 — e 1600 + e" 



w(t) = T J) J mi x . (101) 



/(*) 
J r dtf(t) 

Running the optimization the target functional J% reaches 86.72% of its maximum 



value after 23 iterations. The convergence is presented in figure 6(d) It shows a fast 
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convergence within the first ten iterations. Also shown in the inset of figure 6(b) is the 
increase of the fluence during the iterations. At the end of the iteration the laser fluence 
has reached Eq = 0.7586. 



The optimized laser pulse is shown in figure 6(a) and its spectrum in figure 6(b) 
Obviously, the spectrum shows only the three permitted frequency contributions ujq^ 
W43, and indicated by the vertical lines. 



2.5 




400 
t [a.u.] 

(a) 



800 






r 

': 1 




\l 


■! 




1 


* •'. 


i v 


/A 


j ■ 





C\J 



o 2 

X 

"a 1.5 

CO 



1 

0.5 

0, 



Aid 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 

co [a.u.] 

(b) 



(c) 




Figure 6. (Color online). Control of time-dependent occupation numbers (|0) — > 
|4) — > |3) — * |1 )). (a) Optimized laser field, (b) Spectrum of the optimized pulse 
together with the filter function [( ); scaled by 0.02]. (c) Top panel: Time- 
dependent weight function. Middle panel: Target and the achieved occupations 

\{n\V(t))\ 2 : n = ( ), n = 1 ( ), n = 3 ( ), n = 4 ( ). The 

( ) lines correspond to the target curves |c fe (£)| 2 from equations (|96H100|) . Bottom 

panel: n — 2 ( ), n — 5 (□). (d) Convergence of the algorithm, the ( ) line 

corresponds to the value of the target functional J x , the ( ) line to J, and the 

( ) line to the fluence of the pulse. The vertical line indicates the iteration with 

the maximum yield. 



The system driven by the optimized laser follows the target population, indicated 



Quantum Optimal Control Theory 



33 



by the ( ) lines in figure 6(c), very closely. 

Surprisingly, we observe an occupation of the 5th excited state, indicated by the 
magenta (□) in figure 6(c) , although the resonance frequency U35 is now suppressed by 
the filter. The excitation can only be understood by a mixing of two frequencies, namely 
tuo4 — ^43 = 1.0147 ~ A time-frequency analysis shows that both are present at the 
time the occupation occurs. 

The weight function improves the occupation of the first excited state at the end 
of the pulse. The effect of the weight function on "smoothing" the transition at t — 500 
is rather small compared to runs without the weight function |58j . 

We conclude the analysis with the remark that for the optimizations of time- 
dependent targets we find the choice of the initial guess more important compared 
to the time- independent case. The simplest choice e(t) = turns out to be the best in 
this case, i.e., it produces the highest yields. Other choices end up in similar results 
but slightly smaller yields. The reason is a slowly decaying influence of the initial guess 
during the iteration. 



5. Summary and Outlook 

In this tutorial we have presented the basics of quantum optimal control theory. We 
started with the derivation of the control equations from a suitably formed functional 
and then presented several algorithms to solve these equations. In addition we have 
presented a filter algorithm and the extension to time- dependent targets. Each section 
contained an example to demonstrate the discussed algorithm. 

Currently, we are working on the implementation of several optimal control 
algorithms into the freely available software package OCTOPUS [70j[71]- All algorithms 
in this article are already accessible within OCTOPUS for single-particle calculations 
in real-space. In the future this implementation will allow us to investigate the 
optimal control of multi-electron systems using time- dependent density-functional 
theory (TDDFT)p]. 

Besides further developments in experimental and theoretical capabilities it is 
extremely important for the future of quantum control to work on the interface between 
both "worlds". Namely to extract the most relevant information from the theoretical 
calculations and directly apply them in the experiment. 



Appendix 

Appendix A.l. Closed-loop learning experiments 

In the following we describe the principles of a closed-loop quantum control experiment 
using femtosecond laser pulses. In particular, we will focus on the process of laser pulse 
shaping which provides the connection with the theory presented here. These quantum 
control experiments have become possible due to the improvement of laser pulse 
shaping [3], [72], [73] and the implementation of closed-loop learning (CLL) techniques 
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[B|. Experiments using CLL have delivered highly encouraging results, ranging from the 
control of chemical reactions [U [121 13 El HQl Ell EH] to the control of high-harmonic 
generation [81 [13] . 

Consider a molecule which consists of three parts: A — B — C. The objective is 
to optimize a laser pulse that breaks the bond between A and B, or in other words, 
maximizes the yield of B—C fragments over A—B fragments. The optimization proceeds 
iteratively within a learning loop. The loop starts with using an initial guess for the 
laser pulse. The pulse hits the molecules (in the gas phase) in a reaction chamber. After 
the laser interaction the product is analyzed with a mass-spectrometer. The resulting 
mass spectrum is then fed back into a computer which generates a new pulse shape. 
The prediction of the new pulse shape is usually done by a genetic algorithm [73] . When 
the reaction chamber is loaded with a new sample of the molecule the procedure starts 
again. The loop is continued until the best pulse is found. 

To understand some of the problems in these control experiments we have to take a 
closer look at the process of pulse shaping. The first point to realize is that a femtosecond 
pulse is broad in frequency space, i.e., many frequency components are needed to form 
this pulse (20 fs ~ 1000 cm" 1 ). The idea of pulse shaping is to manipulate the phases 
and the amplitudes of these frequencies which can be achieved in the following way: The 
incoming pulse is targeted on a grating which separates the frequencies in space. Then 
the light beam enters a liquid crystal modulator (LCM or SLM: spatial light modulator) 
which consists of several (typically 128 to 640) small "windows" (pixels). Every pixel 
modifies the amplitude and phase of the incoming light separately, controlled by the 
computer algorithm. The transmitted beam is then transformed back using a second 
grating. The first loop of the experiment has to be started with an initial guess for the 
settings of the LCM. Since the convergence of the genetic algorithm (speed and final 
result) might depend on the initial guess a good choice is very important since the search 
space itself is gigantic: Let us assume a resolution of the amplitude and the phase in 
each LCM pixel of 4 bit which corresponds to 2 4 = 16 different settings. With 128 Pixels 
we have 128^ 2 ' 16 - ) ~ 10 67 different pulses that can be generated. If we fix the amplitude 
setting and use phase-only shaping, we can reduce the number of configurations by a 
factor 1/2. If we give up the idea of free optimization and assume a function which 
describes the settings of the LCM with n parameters [73] we can reduce the search 
space from 128 to n dimensions and therefore n^ 2 ' 16 ) configurations. To determine a 
good parameterization of the pulse, theoretical calculations which take into account the 
experimental constraints are extremely important. 

Appendix A. 2. Two-Level theory 

Two-level systems are applied in many fields of physics from spin models to the area 
of quantum computation [75J. Since they are analytically solvable if the rotating- 
wave approximation (RWA) is used, they often serve as simple but powerful models 
to analyze field-matter interaction problems. The availability of exact solutions is also 
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our motivation to study the system. The idea is to illustrate the OCT algorithm and 
to compare the results to the exact solutions within the RWA. 

First an introduction to the theoretical concepts of two-level systems is given. We 
start by reviewing the RWA and then we derive two recipes for population inversion 
within the RWA (7i~-pulses) and in first-order perturbation theory. We conclude this 



brief review with a controllability analysis in Appendix A. 4 



Appendix A. 3. Two-level system: Basics 

Let our model system consist of two orthonormal states \a) and \b). The state vector 
at any time t is given by 

\*{t)) = ca{t)\a) + c b (t)\b)=: ( ^ V (A.l) 

The time evolution of \^f(t) ) is described by the time- dependent Schrodinger equation 

ift|*(t)) = £(f)|*(t)>, 

with the Hamiltonian in the basis | a ) and | b ) given by 

Hit) ° )- e(t) ( Paa Pab ) , (A.2) 

where we assume p a b = Pb a = A* an< i Paa = Pbb = 0. The time-dependent Schrodinger 
equation then yields the following system of differential equations, 

c a = - iuj a c a (t) + ie(t)pc b (t), (A. 3) 

c b = - iuj b c b {t) + ie(t)pc a {t). (A.4) 

Applying the transformation gk{t) = exp(iukt)ck{t) we obtain 

g a = ie(t) f ie- hJ >» t g b (t), (A.5) 
g b =ie(t)/ie^ a (t), (A.6) 

where u ba = u b - oo a . 

For an arbitrary laser field e(t), the set of differential equations ( 1 A. 51) and ( 1A.6I) is 
only numerically solvable. 

Appendix A. 3.1. Rotating wave approximation (RWA) If the laser field e(t) is chosen 
to be 

e(t) = A sm(ut) = ^{e lvi - e~ [ut ), (A.7) 
we can rewrite the differential equations (1A.5j) and (1A.61) as 
g a (t) = i^e-^*^ (e-* - e~ M ) g b (t) 

= ^R^-i^ ba -u)t _ e -iK +-)t)^( t ) ; (A.8) 
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where we have introduced the expression fi^ = Afi known as the Rabi- frequency. 
Analogously, we get 



9b(t) 



-iA-t „— iA+t 



)»«(*) 



(A.9) 



with A± = (c^v ± z/). If Uf, a and z/ are of similar magnitude (and have the same sign), 
i.e., A_ is small, we may neglect the term 

(A.10) 



e -iA+t _ e -i(u) ba +u)t 



because it oscillates much faster than the other term and cancels out on average. This 
approximation is called the rotating wave approximation (RWA). 



Appendix A. 3. 2. Solutions: Resonant case in RWA Using the RWA, equations (1A.8I) 
and (1A.9j) can be solved analytically [76]. The solution is 

- / ...... iA. 



9a(t) 



g a (0) I cos(f2t/2) Q 



sin(fti/2) 



9b(t) 



-^ 6 (0)sin(f2t/2) 



iA_ 



g b (0) cos(fit/2) 



Mr 



n 



~9a{0) sin(fi*/2) 



jA-t/2 



sin(m/2) 



3 -iA_t/2 



(ah; 



(A.12) 



with 



Q 2 = Q 2 R + A 2 _. 



For the initial state, we choose g a (0) = 1 and g&(0) = 0, i.e., the population is assumed 
to rest in the lower-state at t = 0. Thus, the occupation numbers are given by 



A 2 

\g a {t)\ 2 = cos 2 (fit/2) + — sin 2 (m/2), 
\9b(t)\ 2 = ^|sin 2 (fit/2). 



(A.13) 
(A.14) 



The populations oscillate with the frequency Q/2 which, in the resonant case A_ = 0, 
corresponds to half the Rabi frequency Qr/2. On the other hand, if we fix the 
propagation time to T and look at the final populations for different field amplitudes, 
we also find oscillations of |g a (T)| 2 and |g;,(T)| 2 versus the field amplitude. These are 
called Rabi oscillations as well. 



Appendix A. 3. 3. RWA: Optimal amplitudes If we want to maximize the occupation of 
state | b) at the final time T, it is possible to find an optimal amplitude A by looking at 
the structure of equation (IA.14I) . i.e., if |g;,(T)| 2 reaches its maximum at T it has to be 
zero at 2T: \g b {2T)\ 2 = 0. Therefore, we have 

sin(fi T) = 0, 
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QT = n(2k + 1) with k = 0, 1, ... , 
^ |= <?^_ A >, (A , 5) 

where tells us how many times the maximum occupation has been reached within the 
time interval [0,T]. From equation (IA.15I) follows the optimal amplitude 



„ 1 (2k + l)H 2 A9 /A . 

Note that for a fixed final time T the number k has to be chosen large enough, so that 
the root does not yield an imaginary number for the value of A_. The occupation of 
state | b) for the optimal amplitude in equation (1A.16j) is 



maxW T)P = § = l-p^- 5 . (A.17) 

For a fixed final time T, the result will improve with larger k and therefore with 
increasing field strength [see equation flA.16j) ]. 
The solution in the resonant case (A_ = 0), 

%°A RWA = ^, (A. 18) 

plays an important role in this work since it will be used as a benchmark for optimal 
control solutions. 

The results, equations ( 1A.16I) and ( 1A.18j) . allow us to find a pulse which maximizes 
the occupation in the upper state \b) in an arbitrarily short time. In the limit T — > 
the optimal amplitude diverges as 1/T. However, if T becomes comparable to 1/u the 
RWA is no longer applicable. Besides this limitation, one has to take care that the full 
quantum system can be approximated by a two-level model at all in the case of the 
resulting strong fields, i.e., usually more levels have to be considered. 

Appendix A. 3. 4- Results from perturbation theory In the following we determine an 
optimal laser field for population inversion with the help of perturbation theory. We 
consider the two-level system defined by equations ( 1A.3I) and (1A.4I) . In time-dependent 
perturbation theory [77] the /cth-order coefficients are determined by 

cf ] = -i / ue(t)e i ^'ci fe - 1) , (A. 19) 

c( fc )= -i/ie(t)e- i ^*cf" 1) . (A.20) 

Now consider a functional J = J\ + J2 with 

Ji = \(bMT))\ 2 , (A.21) 

J 2= - A / alt e 2 (t), (A.22) 







where A is a given penalty factor. 

This functional expresses our wish to transfer the occupation, initially in state \a), to 
state \b) at the end of the laser pulse and to minimize the fluence of the laser. 



Quantum Optimal Control Theory 



38 



Setting the variation of this functional to zero yields 
5Jx 



8e{t) 



-2Ae(t). 



(A.23) 



Now we express J\ in first-order perturbation theory using equations flA.19j) 



li I dt' e iUbat ' e(t') 
'o 



(&|*(T))« = c«(T) 
and obtain 

jf) = I|^| 2 f T dt'[ T dt" cos [cu ba (t' - t")} e(t')e(t"). 

The functional derivative of with respect to the laser field is 
5 J? 



|/i| 2 f dt! cos [aj ba (t' - t)] e(t') 
Jo 



Se(t) Jo 
and the optimal field can be determined by 

Xe(t) 



[ dt'\fi\ 2 cos [uj ba {t' -t)]e{t'). 
Jo 



(A.24) 



(A.25) 



(A.26) 



This equation represents an eigenvalue problem which also means that the penalty factor 
A does only yield solutions for certain values, i.e., the eigenvalues of this equation. It can 
be interpreted as the yield per unit of field fluence [78j [79]. Thus, the eigenvector 
corresponding to the largest eigenvalue yields the optimal field (see figure [ATT) . However, 



o 

X 




100 



200 
t[a.u v 



300 



400 



100 



200 
t [a.u.] 



300 



400 



(a) 



(b) 



Figure Al. (Color online), (a) The two largest eigenvectors of equation (|A.26[) . The 

( ) line corresponds to an eigenvalue of 30.8478 while the ( ) line corresponds 

to the eigenvalue 30.8032. The third largest eigenvalue is 1.1126 ■ 10~ 14 . The 
shown eigenvectors are oscillations with uj — u; ba — 0.1568 and have an amplitude 
of A — 0.1273, after rescaling it with equation (|A.27|) . If we propagate the two-level 
system with the field corresponding to the largest eigenvalue we find j[ = 0.6850 

[see (b); ( ) line], showing that the first-order approximation is insufficient in this 

case. The ( ) line corresponds to the ground state occupation. 
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the amplitude of that field is arbitrary. It can be fixed by the maximum value that we 
can achieve for the occupation of state |6), which is one. Therefore, the optimal field is 
given by: 



An optimal field from the two-level system obtained in this way is shown in figure IA1I 
We emphasize that the optimal field obtained in this way is unique and the RWA need 
not be applied. However, as the example calculation shows, the amplitude obtained in 
this way is too small. The final occupation of the target state is 68.50%. 

Appendix A. 4- Example: Complete controllability of two-level-system 

To show that the two-level system is completely controllable we use the theorems of 
Refs. [4"6"t Wl\ , which have been stated in section 12.11 Following these theorems, we have 
to construct the Lie algebra L of the skew-Hermitian operators iH and iHi describing 
our system. In particular, we have to show that the dimension of Lq is N 2 = 2 2 . For the 
construction we will apply the approach of Ref. [17]. The calculation can be simplified 
by rewriting the Hamiltonian of the two-level system ( 1A.2I) as 



where we have used p = p ab = Pba and p aa = pbb = 0. The are the Pauli matrices 
which obey the commutation relation: [o~j, oA = 2x6^^^. At this point we introduce the 
matrix W which will represent the basis of Lq. If its rank is N 2 the system is completely 
controllable. The matrix consists of columns W :> k, i.e., W = (W :j i, W u z, . . . , W.^). The 
columns are calculated by transforming the N x N matrix into a vector with N 2 
components by appending the next column to the end of the previous one. The columns 
are calculated in the following way: 

(i) Set W :! i := Hq to the unperturbed Hamiltonian, here W :> i = (u a , 0, 0, 0Jb)- 

(ii) Set W :> j := Hj to the perturbation/control Hamiltonians, here W :> 2 = (0, ep, ep, Qy. 

(iii) Calculate all non-vanishing commutators of Hq and Hj and add to W if it is 
linearly independent from the existing columns. Here, we have simply W-^ = 
(0,u ba ep, -uj ba ep, 0) f . 

(iv) After appending the new column, all commutators with the previous columns have 
to be evaluated and added (if linearly independent). 

This will result in W :A = {2e 2 p 2 uj ba , 0, 0, ~2e 2 p 2 u ba )\ 

(v) Repeat steps 3 and 4 until the matrix is full, or no new linearly independent columns 
can be found, and then determine the rank of W. 




(A.27) 



iY = -(w a + w i) )/ + - (u a - u b ) o z - e(t)pa. 



(A.28) 
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The two-level system results in the matrix 

/ oo a 2e 2 fi 2 u ba \ 

_ e/i u ba €fi 

efi -uj ba efi 

\ u b -2e 2 /i 2 w 6a J 

If u b ^ —u a the rank of W is N 2 = 4 and therefore the system is completely controllable. 
If uj b = —oj a (i.e. uj ba = — 2u; a ) the rank reduces to iV 2 — 1 = 3 and the system is not 
completely controllable according to the definition. In other words, the difference of 
being controllable or not depends on whether Hq has a non-zero trace or not. The 
disturbing fact is that a traceless Ho can be changed to a Hamiltonian with non-zero 
trace by shifting the energy levels. However, this shift is not of physical significance. 
The dilemma is resolved by realizing that in the case of a traceless Hq we (only) loose 
control over the phase of the state, which is not relevant in the case of population 
inversion [47] . 

Acknowledgments 

We thank Elham Khosravi for reviewing the manuscript. This work was supported, 
in part, by the Deutsche Forschungsgemeinschaft through the SFB450, and the 
NANOQUANTA Network of Excellence of the European Union. 

References 

[1] H. Rabitz, R. de Vivie-Riedle, M. Motzkus, and K. Kompa, Science 288, 824 (2000). 
[2] W. S. Warren, H. Rabitz, and M. Dahleh, Science 259, 1581 (1993). 

[3] A. M. Weiner, D. E. Leaird, J. S. Patel, and J. R. Wullert, IEEE J. Quant. Electron. 28, 908 
(1992). 

[4] A. Assion, T. Baumert, M. Bergt, T. Brixner, B. Kiefer, V. Seyfried, M. Strehle, and G. Gerber, 

Science 282, 919 (1998). 
[5] R. J. Levis, G. M. Menkir, and H. Rabitz, Science 292, 709 (2001). 
[6] R. S. Judson and H. Rabitz, Phys. Rev. Lett. 68, 1500 (1992). 

[7] C. J. Bardeen, V. V. Yakovlev, K. R. Wilson, S. D. Carpenter, P. M. Weber, and W. S. Warren, 

Chem. Phys. Lett. 280, 151 (1997). 
[8] R. Bartels, S. Backus, E. Zeek, L. Misoguti, G. Vovin, I. P. Christov, M. M. Murmane, and H. C. 

Kapteyn, Nature 406, 164 (2000). 
[9] T. Brixner, N. H. Damrauer, P. Niklaus, and G. Gerber, Nature 414, 57 (2001). 
[10] J. L. Herek, D. Z. W. Wohlleben, R. J. Cogdell, and M. Motzkus, Nature 417, 533 (2002). 
[11] T. Brixner, N. H. Damrauer, B. Kiefer, and G. Gerber, J. Chem. Phys. 118, 3692 (2003). 
[12] C. Daniel, J. Full, L. Gonzalez, C. Lupulescu, J. Manz, A. Merli, S. Vajda, and L. Woste, Science 
299, 536 (2003). 

[13] T. Pfeifer, D. Walter, C. Winterfeldt, C. Spielmann, and G. Gerber, Appl. Phys. B 80, 277 (2005). 
[14] G. Vogt, G. Krampert, P. Niklaus, P. Nuernberger, and G. Gerber, Phys. Rev. Lett. 94, 068305 
(2005). 

[15] L. Polachek, D. Oron, and Y. Silberberg, Opt. Lett. 31, 5 (2006). 

[16] M. Plewicki, S. M. Weber, F. Weise, and A. Lindinger, Appl. Phys. B. (2006), submitted. 
[17] M. Plewicki, F. Weise, S. M. Weber, and A. Lindinger, Appl. Opt. (2006), submitted. 



Quantum Optimal Control Theory 



41 



[18] G. M. Huang, T. J. Tarn, and J. W. Clark, J. Math. Phys. 24, 2608 (1983). 

[19] B. Schafer-Bung, R. Mitric, V. Bonacic-Koutecky, A. Bartelt, C. Lupulcscu, A. Lindinger, S. Vajda, 

S. M. Weber, and L. Woste, J. Phys. Chem. A 108, 4175 (2004). 

[20] W. Jakubetz, B. Just, J. Manz, and H. J. Schreier, J. Phys. Chem. 94, 2294 (1990). 

[21] D. J. Tannor and S. A. Rice, J. Chem. Phys. 83, 5013 (1986). 

[22] D. J. Tannor, R. Kosloff, and S. A. Rice, J. Chem. Phys. 85, 5805 (1986). 

[23] J. P. R. Mitric, M. Hartmann and V. Bonacic-Koutecky, Eur. Phys. J. D 24, 177 (2003). 

[24] P. Brumcr and M. Shapiro, Chem. Phys. Lett. 126, 541 (1986). 

[25] U. Gaubatz, P. Rudecki, S. Schiemann, and K. Bcrgmann, J. Chem. Phys. 92, 5363 (1990). 

[26] G. W. Coulston and K. Bcrgmann, J. Chem. Phys. 96, 3467 (1992). 

[27] K. Bcrgmann, H. Thcuer, and B. W. Shore, Rev. Mod. Phys. 70, 1003 (1998). 

[28] X. Chu and S. I. Chu, Phys. Rev. A 64, 021403 (2001). 

[29] R. Kosloff, S. A. Rice, P. Gaspard, S. Tersigni, and D. J. Tannor, Chem. Phys. 139, 201 (1989). 

[30] A. P. Peirce, M. A. Dahleh, and H. Rabitz, Phys. Rev. A 37, 4950 (1988). 

[31] D. J. Tannor, V. Kazakov, and V. Orlov, Time- Dependent Quantum Molecular Dynamics, chapter 
Control of photochemical branching: Novel procedures for finding optimal pulses and global 
upper bounds, p. 347, Plenum Press, New York, 1992. 

[32] W. Zhu, J. Botina, and H. Rabitz, J. Chem. Phys. 108, 1953 (1998). 

[33] W. Zhu and H. Rabitz, J. Chem. Phys. 109, 385 (1998). 

[34] Y. Maday and G. Turinici, J. Chem. Phys. 118, 8191 (2003). 

[35] A. Bartana, R. Kosloff, and D. J. Tannor, J. Chem. Phys. 99, 196 (1993). 

[36] A. Bartana, R. Kosloff, and D. J. Tannor, J. Chem. Phys. 106, 1435 (1997). 

[37] Y. Ohtsuki, W. Zhu, and H. Rabitz, J. Chem. Phys. 110, 9825 (1999). 

[38] Y. Ohtsuki, K. Nakagami, Y. Fujimura, W. Zhu, and H. Rabitz., J. Chem. Phys. 114, 8867 (2001). 

[39] S. Lloyd, Phys. Rev. A 62, 022108 (2000). 

[40] N. Weaver, J. Math. Phys. 41, 5262 (2000). 

[41] N. Khancja, R. Brockett, and S. J. Glaser, Phys. Rev. A 63, 032308 (2001). 

[42] R. Wu, C. Li, and Y. Wang, Physics Letters A 295, 30 (2002). 

[43] H. Jirari and W. P6tz, Phys. Rev. A 72, 013409 (2005). 

[44] S. Beyvers, Y. Ohtsuki, and P. Saalfrank, J. Chem. Phys. 124, 234706 (2006). 

[45] T. C. Weinacht and P. H. Bucksbaum, J. Opt. B: Quantum Semiclass. Opt. 4, R35 (2001). 

[46] V. Ramakrishna, M. V. Salapaka, M. Dahleh, H. Rabitz, and A. Peirce, Phys. Rev. A 51, 960 
(1995). 

[47] S. G. Schirmcr, H. Fu, and A. I. Solomon, Phys. Rev. A 63, 063410 (2001). 

[48] K. Sundcrmann and R. de Vivie-Ricdle, J. Chem. Phys. 110, 1896 (1999). 

[49] G. Turinici and H. Rabitz, Chem. Phys. 267, 1 (2001). 

[50] M. Holthaus and B. Just, Phys. Rev. A 49, 1950 (1994). 

[51] J. Werschnik and E. K. U. Gross, J. Opt. B: Quantum and Semiclass. Opt. 7, S300 (2005). 

[52] M. Grifoni and P. Hanggi, Physics Reports 304, 229 (1998). 

[53] J. A. Fleck, J. R. Morris, and M. D. Fcit, Appl. Phys. A 10, 129 (1976). 

[54] M. Frigo and S. G. Johnson, FFTW: An adaptive software architecture for the FFT, in Proc. 

1998 IEEE Intl. Conf. Acoustics Speech and Signal Processing, volume 3, pp. 1381-1384, Seattle, 

WA , USA, 1998, IEEE. 

[55] P. Gross, D. Neuhauser, and H. Rabitz, J. Chem. Phys. 96, 2834 (1992). 

[56] I. R. Sola, J. Santamaria, and D. J. Tannor, J. Phys. Chem. A 102, 4301 (1998). 

[57] T. Hornung, M. Motzkus, and R. de Vivie-Riedle, J. Chem. Phys 115, 3105 (2001). 

[58] J. Werschnik, Quantum Optimal Control Theory: Filter Techniques, Time- Dependent Targets, 

and Time- Dependent Density- Functional Theory, Cuvillier Verlag, Gottingen, 2005. 

[59] S. Shi and H. Rabitz, J. Chem. Phys. 92, 364 (1990). 

[60] I. Grigorenko, M. Garcia, and K. H. Benncmann, Phys. Rev. Lett. 89, 233003 (2002). 

[61] W. Zhu and H. Rabitz, J. Chem. Phys. 119, 3619 (2003). 



Quantum Optimal Control Theory 



42 



2] M. Sugawara, J. Chem. Phys. 118, 6784 (2003). 

3] Y. Ohtsuki, G. Turinici, and H. Rabitz, J. Chem. Phys. 120, 5509 (2004). 
4] I. Serban, J. Werschnik, and E. K. U. Gross, Phys. Rev. A 71, 053810 (2005). 
5] A. Kaiser and V. May, J. Chem. Phys. 121, 2528 (2004). 

3] I. Serban, Optimal control of time-dependent targets, Master's thesis, Freie Universitat Berlin, 
2004. 

[67] J. Werschnik and E. K. U. Gross, Spie Proceedings (2006). 

[68] J. J. Sakurai, Modern Quantum Mechanics, Addison Wesley, Reading, Massachusetts, revised 
edition, 1994. 

[69] A. Kaiser and V. May, Chem. Phys. Lett. 405, 339 (2005). 

[70] M. A. L. Marques, A. Castro, G. F. Bertsch, and A. Rubio, Comput. Phys. Commun. 151, 60 
(2003). 

[71] http://www.tddft.org/octopus . 

[72] T. Brixncr, G. Krampert, T. Pfeifer, R. Selle, G. Gerber, M. Wollenhaupt, O. Graefe, C. Horn, 

D. Liese, and T. Baumcrt, Phys. Rev. Lett. 92, 208301 (2004). 
[73] M. Y. Shverdin, D. R. Walker, D. D. Yavuz, G. Y. Yin, and S. E. Harris, Phys. Rev. Lett. 94, 

033904 (2005). 

[74] D. Zcidlcr, S. Frcy, K. L. Kompa, and M. Motzkus, Phys. Rev. A 64, 023420 (2001). 
[75] M. A. Nielsen and I. L. Chuang, Quantum Computation, Cambridge University Press, Cambridge, 
2004. 

[76] M. O. Scully and M. S. Zubairy, Quantum Optics, Cambridge University Press, Cambridge, U.K., 
1997. 

[77] B. H. Bransden and C. J. Joachain, Physics of atoms and molecules, Longman, New York, 1983. 
[78] J. L. Krause, R. M. Whitnell, K. R. Wilson, Y. Yan, and S. Mukamcl, J. Chem. Phys. 99, 6562 
(1993). 

[79] J. L. Krause, M. Messina, K. R. Wilson, and Y. Yan, J. Phys. Chem. 99, 13736 (1995). 



